Standard library vectorised functions

for (i in 1:nrow(df)) {
  # New column based on 2 others
  if (df$x[i] > 5 && df$y[i] < 0) {
    df$y[i] <- 5
  } else {
    df$y[i] <- 0
  }
  
  # Replace NAs with error code
  if (is.na(df$z[i])) {
    df$z[i] <- 9999 
  }
  
  # String concatenate columns
  df$xy[i] <- paste(df$x[i], df$y[i], sep="_")
}
  
# Distance between every row
dists <- matrix(nrow=nrow(df), ncol=nrow(df))
for (i in 1:nrow(df)) {
  for (j in 1:nrow(df)) {
    if (j != i) {
      dists[i, j] <- sqrt((df$x[i] - df$x[j])**2 + (df$y[i] - df$y[j])**2)
    }
  }
}

# New column based on 2 others
df$y <- ifelse(df$x > 5 & df$y < 0, 5, 0)





# Replace NAs with error code
df$z[is.na(df$z)] <- 9999



# Concatenate columns
df$xy <- paste(df$x, df$y, sep="_")


# Distance between every row
dist(df[, c('x', 'y')])

Conclusions

In general, if you’re using a for loop (or a sapply variant), then your code could be sped up by using a vectorised function - Me (7 slides ago)

• This wasn’t fully correct, as vectorised functions can have for loops under the hood and will thus still be slow
• The difference between a vectorised function built using Vectorize or an inbuilt function like basename is that the latter will have a for loop, but it will be written in C/C++ rather than R

Basic DataFrame operations

• Fortunately working with data.frames and the tidyverse core verbs pushes you towards using vectorised functions
• group_by() |> summarise() is both quicker and more legible than manually looping over the groups and combining the results
• filter(f(x)) assumes that f() is vectorised and returns a boolean TRUE/FALSE for every row
• mutate(newcol=f(oldcol)) assumes f() is vectorised and returns a value per row

# Non-vectorised version didn't Error, 
# but gave an unexpected result
data.frame(path=paths) |>
  mutate(path_clean1 = example_1(path),
         path_clean2 = example_3(path)) |>
  head()

                           path path_clean1 path_clean2
1         some/path/abc/001.txt xyz_002.txt abc_001.txt
2 another/directory/xyz/002.txt xyz_002.txt xyz_002.txt
3         some/path/abc/001.txt xyz_002.txt abc_001.txt
4 another/directory/xyz/002.txt xyz_002.txt xyz_002.txt
5         some/path/abc/001.txt xyz_002.txt abc_001.txt
6 another/directory/xyz/002.txt xyz_002.txt xyz_002.txt

# This function isn't vectorised due to the if/else statement
# Solution: Use ifelse() instead
replace_both_NA_9999 <- function(x, y) {
  if (is.na(x) && is.na(y)) {
    return(9999)
  } else {
    return(0)
  }
}

data.frame(a=c(5, 3, NA, 2, NA), 
           b=c(NA, 2, NA, 1, 9)) |>
  mutate(c = replace_both_NA_9999(a, b))

Error in `mutate()`:
ℹ In argument: `c = replace_both_NA_9999(a, b)`.
Caused by error in `is.na(x) && is.na(y)`:
! 'length = 5' in coercion to 'logical(1)'

Joining

• Linking 2 datasets together using the join family of functions is an integral part of data analysis
• However, join functions are highly efficient functions and can be useful in a number of siutations, even when we don’t have 2 separate datasets
• inner_join links two dataframes together based on a column in common, with the number of rows equal to the number of rows in the ‘left’ table that have a matching row in the ‘right’ table

joined <- df_1 |> inner_join(df_2, by="group")

Example usage: inner join instead of ifelse

• Can think of inner_join as being able to both filter and mutate new columns
• Example: apply different per-group scaling factor to 300,000 measurements from 3 groups
• On one joining column join isn’t much quicker, but it’s far more legible and scales well to both having more groups in the joining column, and additional joining columns

f_join <- function() {
  df |> inner_join(scales, by="group")
}

f_ifelse <- function() {
  df |>
    mutate(scale = ifelse(group == 'a', 2,
                          ifelse(group == 'b', 7.8,
                                 ifelse(group == 'c', 9, NA))))
  
}

res <- microbenchmark(f_join(), f_ifelse(), times=10)
summary(res)[c("expr", "median")]

        expr   median
1   f_join() 22.45482
2 f_ifelse() 27.31786

Interval joins

• Joins aren’t limited to joining on equal values, can also join on intervals or closest value
• Example: Have measurements from every day in 2020, but want to limit analysis to 5 specific weeks

joined <- df_interval |>
  inner_join(weeks, 
             by=join_by(time >= week_start, time < week_end))

Benchmark

• On only 366 rows with 5 groups it is 10x as fast, will scale better, and is more understandable

f_intervaljoin <- function() {
  df_interval |>
    inner_join(weeks, by=join_by(time >= week_start, time < week_end))
}

f_ifelse <- function() {
  df_interval |>
    mutate(week_group = ifelse(time >= as_date("2020-02-14") & time < as_date("2020-02-21"),
                               'a',
                               ifelse(time >= as_date("2020-03-17") & time < as_date("2020-03-24"),
                                      'b',
                                      ifelse(time >= as_date("2020-05-08") & time < as_date("2020-05-15"),
                                             'c',
                                             ifelse(time >= as_date("2020-09-20") & time < as_date("2020-09-27"),
                                                    'd',
                                                    ifelse(time >= as_date("2020-11-13") & time < as_date("2020-11-20"),
                                                           'e', 
                                                           NA)))))) |>
    filter(!is.na(week_group))
}

res <- microbenchmark(f_intervaljoin(), f_ifelse(), times=10)
summary(res)[c("expr", "median")]

              expr    median
1 f_intervaljoin()  2.254611
2       f_ifelse() 27.822321

Example dataset

• What if we’re using fast functions but still experiencing slow performance due to dataset’s size?
• Example dataset: Company House data containing 5 million rows (440MB archive download, extracts to 2.4GB) of all companies incorporated in the UK since 1856
• Using first million rows as an example

df <- read_csv("BasicCompanyDataAsOneFile-2023-05-01.csv", n_max=1e6, show_col_types=FALSE)
df$IncorporationDate <- as_date(df$IncorporationDate, format="%d/%m/%Y")
dim(df)

[1] 1000000      55

df |> 
  select(CompanyName, RegAddress.PostTown, IncorporationDate, SICCode.SicText_1) |>
  head()

# A tibble: 6 × 4
  CompanyName            RegAddress.PostTown IncorporationDate SICCode.SicText_1
  <chr>                  <chr>               <date>            <chr>            
1 ! HEAL UR TECH LTD     GUILDFORD           2022-10-12        33140 - Repair o…
2 ! LTD                  LEEDS               2012-09-11        99999 - Dormant …
3 !? LTD                 ROMILEY             2018-06-05        47710 - Retail s…
4 !BIG IMPACT GRAPHICS … LONDON              2018-12-28        18129 - Printing…
5 !GOBERUB LTD           BISHOP'S STORTFORD  2021-05-17        62020 - Informat…
6 !NFOGENIE LTD          LONDON              2021-07-21        58290 - Other so…

Question 3: Classifications

• Companies can be assigned with up to 4 classifications from a list of 1,042 options
• Do classifications tend to cluster together? I.e. is the average number of classifications a company has related to the first classification?
• Slightly tenuous example but wanted to demonstrate pivoting + joining!
• Only want to look at classifications that are used by at least 10 companies (inner_join to filter)
• Multiple classifications are stored in 4 wide columns that are NA when unused - easier to count the number of non-null column entries in long format

# A tibble: 6 × 5
  CompanyName              SICCode.SicText_1 SICCode.SicText_2 SICCode.SicText_3
  <chr>                    <chr>             <chr>             <chr>            
1 ! HEAL UR TECH LTD       33140 - Repair o… 47421 - Retail s… <NA>             
2 ! LTD                    99999 - Dormant … <NA>              <NA>             
3 !? LTD                   47710 - Retail s… <NA>              <NA>             
4 !BIG IMPACT GRAPHICS LI… 18129 - Printing… 59112 - Video pr… 63120 - Web port…
5 !GOBERUB LTD             62020 - Informat… 70229 - Manageme… 79110 - Travel a…
6 !NFOGENIE LTD            58290 - Other so… <NA>              <NA>             
# ℹ 1 more variable: SICCode.SicText_4 <chr>

Question 3: Classifications (code)

# 755 rows containing the SIC codes that at least 10 companies have
# Only 1 column, SICCode.SicText_1
sic_10_companies <- df |> 
                count(SICCode.SicText_1) |>
                filter(n >= 10) |>
                select(SICCode.SicText_1)

df |>
  # Could do a filter to restrict to these 10 companies, but it's actually quicker to use an inner join
  inner_join(sic_10_companies, by="SICCode.SicText_1") |>
  select(CompanyNumber, SICCode.SicText_1, SICCode.SicText_2, SICCode.SicText_3, SICCode.SicText_4) |> 
  mutate(first_classification = SICCode.SicText_1) |>
  # Pivoting to make it easier to count how many non-NULL classifications each company has
  pivot_longer(c(SICCode.SicText_1, SICCode.SicText_2, SICCode.SicText_3, SICCode.SicText_4)) |>
  filter(!is.na(value)) |>
  # Count how many classifications each company has
  count(CompanyNumber, first_classification) |>
  # Calculate the average number per the first classification
  group_by(first_classification) |>
  summarise(mean_classifications = mean(n, na.rm=T)) |>
  arrange(desc(mean_classifications)) |>
  head()

# A tibble: 6 × 2
  first_classification                                      mean_classifications
  <chr>                                                                    <dbl>
1 7210 - Hardware consultancy                                               3   
2 07100 - Mining of iron ores                                               2.48
3 10611 - Grain milling                                                     2.43
4 18110 - Printing of newspapers                                            2.43
5 14131 - Manufacture of other men's outerwear                              2.40
6 01280 - Growing of spices, aromatic, drug and pharmaceut…                 2.39

Introduction

• data.table is an alternative to data.frame/tibble that is optimised for speed and low memory usage
• The trade-off is that its API is a bit/lot less user friendly

library(data.table)
dt <- fread("BasicCompanyDataAsOneFile-2023-05-01.csv", nrows=1e6)         # fread is the equivalent of read.csv
dt[, IncorporationDate := as_date(IncorporationDate, format="%d/%m/%Y") ]  # Creates a new column by *reference*
dim(dt)

[1] 1000000      55

# Display rows 1-5 and the specified columns
dt[1:5, .(CompanyName, RegAddress.PostTown, IncorporationDate, SICCode.SicText_1)]

                    CompanyName RegAddress.PostTown IncorporationDate
1:           ! HEAL UR TECH LTD           GUILDFORD        2022-10-12
2:                        ! LTD               LEEDS        2012-09-11
3:                       !? LTD             ROMILEY        2018-06-05
4: !BIG IMPACT GRAPHICS LIMITED              LONDON        2018-12-28
5:                 !GOBERUB LTD  BISHOP'S STORTFORD        2021-05-17
                                       SICCode.SicText_1
1:                33140 - Repair of electrical equipment
2:                               99999 - Dormant Company
3: 47710 - Retail sale of clothing in specialised stores
4:                               18129 - Printing n.e.c.
5: 62020 - Information technology consultancy activities

Number of classifications

• Joins are less intuitive. x[y] is equal to left_join(y, x), NOT inner_join(x, y)
• melt is equivalent to pivot_longer and IMO less intuitive
• Intermdiate variables everywhere!

sic_10_companies_dt <- dt[, .N, by=.(SICCode.SicText_1)][ N >= 10, .(SICCode.SicText_1) ]
dt_companies_wide <- dt[ sic_10_companies_dt,  # This is a join!
                         .(CompanyNumber, 
                           first_classification = SICCode.SicText_1,
                           SICCode.SicText_1,
                           SICCode.SicText_2,
                           SICCode.SicText_3,
                           SICCode.SicText_4),
                          on=.(SICCode.SicText_1)]
dt_companies_long <- melt(dt_companies_wide, id.vars=c('CompanyNumber', 'first_classification'))
dt_companies_mean <- dt_companies_long[ value != '',  # Removes the unused SIC columns
                                        .N, 
                                        by=.(CompanyNumber, first_classification)][, 
                                                                                   .(mean_classifications = mean(N, na.rm=T)), 
                                                                                   by=.(first_classification)]
head(dt_companies_mean[ order(mean_classifications, decreasing = TRUE)])

                                                 first_classification
1:                                        7210 - Hardware consultancy
2:                                        07100 - Mining of iron ores
3:                                     18110 - Printing of newspapers
4:                                              10611 - Grain milling
5:                       14131 - Manufacture of other men's outerwear
6: 01280 - Growing of spices, aromatic, drug and pharmaceutical crops
   mean_classifications
1:             3.000000
2:             2.478261
3:             2.428571
4:             2.428571
5:             2.401606
6:             2.392157

tidytable: introduction

library(tidytable)
# Here we explicitly create tidytable from a regular data.frame
# But passing a regular data.frame or data.table into any tidytable function
# will implicitly change it to be a tidytable object
dtt <- as_tidytable(df)

dtt |> 
    count(SICCode.SicText_1) |>
    filter(n >= 10) |>
    select(SICCode.SicText_1) 

# A tidytable: 755 × 1
   SICCode.SicText_1                                                       
   <chr>                                                                   
 1 01110 - Growing of cereals (except rice), leguminous crops and oil seeds
 2 01120 - Growing of rice                                                 
 3 01130 - Growing of vegetables and melons, roots and tubers              
 4 01160 - Growing of fibre crops                                          
 5 01190 - Growing of other non-perennial crops                            
 6 01210 - Growing of grapes                                               
 7 01220 - Growing of tropical and subtropical fruits                      
 8 01240 - Growing of pome fruits and stone fruits                         
 9 01250 - Growing of other tree and bush fruits and nuts                  
10 01270 - Growing of beverage crops                                       
# ℹ 745 more rows

dtplyr: introduction

library(dtplyr)

# dtplyr operates on `lazy data.tables` which are only created by this function
dtp <- lazy_dt(df)

dtp |> 
    count(SICCode.SicText_1) |>
    filter(n >= 10) |>
    select(SICCode.SicText_1) 

Source: local data table [755 x 1]
Call:   `_DT1`[, .(n = .N), keyby = .(SICCode.SicText_1)][n >= 10, .(SICCode.SicText_1)]

  SICCode.SicText_1                                                       
  <chr>                                                                   
1 01110 - Growing of cereals (except rice), leguminous crops and oil seeds
2 01120 - Growing of rice                                                 
3 01130 - Growing of vegetables and melons, roots and tubers              
4 01160 - Growing of fibre crops                                          
5 01190 - Growing of other non-perennial crops                            
6 01210 - Growing of grapes                                               
# ℹ 749 more rows

# Use as.data.table()/as.data.frame()/as_tibble() to access results

dtplyr: usage

dtp |> 
    count(SICCode.SicText_1) |>
    filter(n >= 10) |>
    select(SICCode.SicText_1) |>
    show_query()

`_DT1`[, .(n = .N), keyby = .(SICCode.SicText_1)][n >= 10, .(SICCode.SicText_1)]

dtp |> 
    count(SICCode.SicText_1) |>
    filter(n >= 10) |>
    select(SICCode.SicText_1) |>
    collect() |> 
    head()

# A tibble: 6 × 1
  SICCode.SicText_1                                                       
  <chr>                                                                   
1 01110 - Growing of cereals (except rice), leguminous crops and oil seeds
2 01120 - Growing of rice                                                 
3 01130 - Growing of vegetables and melons, roots and tubers              
4 01160 - Growing of fibre crops                                          
5 01190 - Growing of other non-perennial crops                            
6 01210 - Growing of grapes                                               

dtplyr: chaining queries

• dtplyr queries that haven’t been collect() can be used in joins

# NB: this returns a datatable QUERY, not a dataset itself
sic_10_companies_dtp <- dtp |> 
    count(SICCode.SicText_1) |>
    filter(n >= 10) |>
    select(SICCode.SicText_1) 

# Can join that query into the middle of another query to return another query
results_dtp <- dtp |>
  inner_join(sic_10_companies_dtp, by="SICCode.SicText_1") |>
  select(CompanyNumber, SICCode.SicText_1, SICCode.SicText_2, SICCode.SicText_3, SICCode.SicText_4) |> 
  mutate(first_classification = SICCode.SicText_1) |>
  pivot_longer(c(SICCode.SicText_1, SICCode.SicText_2, SICCode.SicText_3, SICCode.SicText_4)) |>
  filter(!is.na(value)) |>
  count(CompanyNumber, first_classification) |>
  group_by(first_classification) |>
  summarise(mean_classifications = mean(n, na.rm=T)) |>
  arrange(desc(mean_classifications))

# Finally execute the full query
results_dtp |>
  collect() |>
  head()

# A tibble: 6 × 2
  first_classification                                      mean_classifications
  <chr>                                                                    <dbl>
1 7210 - Hardware consultancy                                               3   
2 07100 - Mining of iron ores                                               2.48
3 10611 - Grain milling                                                     2.43
4 18110 - Printing of newspapers                                            2.43
5 14131 - Manufacture of other men's outerwear                              2.40
6 01280 - Growing of spices, aromatic, drug and pharmaceut…                 2.39

Introduction

• All these options require reading the full dataset into memory, not viable if we have larger than memory data
• Embedded relational databases are stored on disk and only read into memory as needed

SQLite usage

# tbl requires a connection and a table name to read data from
tbl(con_sql, "data_1e6") |> 
  filter(RegAddress.PostTown == 'YORK') |> 
  mutate(postcode = word(RegAddress.PostCode, 1)) |>
  count(postcode) |>
  arrange(desc(n)) |>
  head(5) |>
  show_query()

<SQL>
SELECT `postcode`, COUNT(*) AS `n`
FROM (
  SELECT *, word(`RegAddress.PostCode`, 1.0) AS `postcode`
  FROM `data_1e6`
  WHERE (`RegAddress.PostTown` = 'YORK')
)
GROUP BY `postcode`
ORDER BY `n` DESC
LIMIT 5

tbl(con_sql, "data_1e6") |> 
  filter(RegAddress.PostTown == 'YORK') |> 
  mutate(postcode = substr(RegAddress.PostCode, 1, 4)) |>
  count(postcode) |>
  arrange(desc(n)) |>
  head(5) |>
  collect()

# A tibble: 5 × 2
  postcode     n
  <chr>    <int>
1 "YO30"     486
2 "YO19"     325
3 "YO26"     271
4 "YO1 "     241
5 "YO31"     180

duckdb: usage

tbl(con_dd, "data_1e6") |> 
  filter(RegAddress.PostTown == 'YORK') |> 
  mutate(postcode = substr(RegAddress.PostCode, 1, 4)) |>
  count(postcode) |>
  arrange(desc(n)) |>
  head(50) |>
  show_query()

<SQL>
SELECT postcode, COUNT(*) AS n
FROM (
  SELECT *, SUBSTR("RegAddress.PostCode", 1, 4) AS postcode
  FROM data_1e6
  WHERE ("RegAddress.PostTown" = 'YORK')
) q01
GROUP BY postcode
ORDER BY n DESC
LIMIT 50

tbl(con_dd, "data_1e6") |> 
  filter(RegAddress.PostTown == 'YORK') |> 
  mutate(postcode = substr(RegAddress.PostCode, 1, 4)) |>
  count(postcode) |>
  arrange(desc(n)) |>
  head(50) |>
  collect()

# A tibble: 23 × 2
   postcode     n
   <chr>    <dbl>
 1 "YO30"     486
 2 "YO19"     325
 3 "YO26"     271
 4 "YO1 "     241
 5 "YO31"     180
 6 "YO24"     178
 7 "YO10"     172
 8 "YO32"     169
 9 "YO23"     161
10 "YO42"     142
# ℹ 13 more rows

Overall benchmark

Introduction

• Sometimes for loops are necessary:
  • No inbuilt vectorised solution
  • Recurrent algorithm such as stepping through time or space
  • Performant critical code and need more specialised data structures
• Rcpp, which combines R and C++, to the rescue!

Basic example

• Can use Rcpp::cppFunction() to write a C++ function as a string or Rcpp::sourceCpp() if it’s in a separate file
• Both methods do the same:
  • Compile C++ code
  • Create an R function that calls it
• C++ has many differences with R, having to assign every variable a type is probably the most notable

library(Rcpp)
cppFunction("
double sumRcpp(NumericVector x) {
  // The function definition syntax is:
  // RETURN_TYPE functionName(INPUT_TYPE input, ...)
  int n = x.size();  // R objects have their own type (NumericVector) with useful attributes
  double total = 0;  // Need to instantiate variables before use
  for (int i = 0; i < n; i++) {  // C++ indexes start at 0
    total += x[i];
  }
  return total;      // Need to explicitly return values
}")
# The sumRcpp function is instantly available within R
sumRcpp(c(1, 2, 3))

[1] 6

Syntatic sugar - data structures

• When calling an Rcpp function, the inputs are automatically converted from their R type into the specified C++ type
• NumericVector is a special Rcpp data structure that represents a vector of floats, so will accept both c(1.2, 2.4, 3.6), and c(1, 2, 3), but not c('a', 'b', 'c')
• IntegerVector coerces floats into integers
• CharacterVector also exists for strings, and NumericMatrix for 2D structures

sum benchmarks

sumR <- function(x) {
  total <- 0
  for (i in seq(length(x))) {
    total <- total + x[i]
  }
  total
}

Real world example: Kalman Filter

kf_r <- function(y, m, Q=0.5, H=0.5) {
  n <- length(y)
  alpha <- array(NA, dim=c(n+1, m))
  P <- array(NA, dim=c(m, m, n+1))
  alpha[1] <- 0  # Initialise with zero mean and high variance
  P[, , 1] <- 1e3
  Z <- array(1, dim=c(n, m)) 
  for (i in 1:n) {
    P_updated <- P[, , i] + Q
    # Calculate kalman gain
    K <- P_updated %*% t(Z[i, ]) %*% solve(Z[i, ] %*% P_updated %*% t(Z[i, ]) + H)
    # Update state and covariance
    alpha[i+1, ] <- alpha[i, ] + K %*% (y[i] - Z[i, ] %*% alpha[i, ])
    P[, , i+1] <- (diag(m) - K %*% Z[i, ]) %*% P_updated
  }
  list(alpha=alpha, P=P)
}

Kalman Filter - Rcpp implementation

cppFunction("
NumericVector kf_rcpp(arma::vec y, int m, float Q=0.5, float H=0.5) {
  int n = y.n_rows;
  
  arma::mat alpha(n+1, m, arma::fill::none);
  arma::cube P(m, m, n+1, arma::fill::none);
  arma::mat Z(n, m, arma::fill::ones);
  
  // Initialise with zero mean and high variance
  alpha.row(0).fill(0);
  P.slice(0).diag().fill(1000);
  
  // Run filter
  arma::mat P_updated(m, m);
  arma::mat K(m, m);
  for (int i=0; i<n; i++) {
    P_updated = P.slice(i) + Q;
    // Calculate kalman gain:
    K = P_updated * Z.row(i).t() * arma::inv(Z.row(i) * P_updated * Z.row(i).t() + H);
    // Update state and covariance
    alpha.row(i+1) = alpha.row(i) + K * (y[i] - Z.row(i) * alpha.row(i));
    P.slice(i+1) = (arma::eye(m, m) - K * Z.row(i)) * P_updated;
  }
  
  return wrap(alpha);  // This is crucial, converts the Armadillo matrix into an R NumericVector
}", depends="RcppArmadillo")

Parallelisation / Viking

• For iterative jobs that don’t fit within tabular data can run parallelised for loops, e.g.:
  • Fitting models
  • Network downloads/uploads
  • File processing
• For ‘small’ jobs can run locally using parallel::mclapply (Linux), doParallel and foreach (Windows), or furrr (all OS, Tidyverse, combines future and purrr)

Thoughts

• I tend to use all of these strategies with varying frequency:
  • Using inbuilt vectorised functions - daily
  • data.table in some form - weekly
  • Parellel jobs (Viking or local) - monthly
  • Rcpp - couple of times a year
• As well as speeding up interactive work, faster code is especially useful when developing packages
• As a result of researching for this talk, I’m going to use duckdb in future for some large datasets

Resources

• Vectorisation: Chapter in Advanced R
• Joining: Tidyverse docs, interactive visual join viewer (nb: ‘outer’ joins are called ‘full’ joins in R)
• data.table: vignette, syntax comparison with tidyverse
• tidytable vs dtplyr: benchmarking between tidytable, data.table, dtplyr, and tidyverse and pandas (nb: from tidytable author)
• SQLite: tutorial
• duckdb: official docs
• Rcpp: chapter in Advanced R, Rcpp book (thorough), Rcpp for everyone book (accessible)
• Parallelisation: chapter in R Programming for Data Science ebook
• Viking: wiki